home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / ibeta_pdf.pro < prev    next >
Text File  |  1997-07-08  |  2KB  |  62 lines
  1. ;$Id: ibeta_pdf.pro,v 1.3 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1994-1997, Research Systems, Inc.  All rights reserved.
  4. ;       Unauthorized reproduction prohibited.
  5. ;+
  6. ; NAME:
  7. ;       IBETA_PDF
  8. ;
  9. ; PURPOSE:
  10. ;       This function computes the incomplete beta function.
  11. ;       It is called by the probability density functions in 
  12. ;       this directory. See the function IBETA() in the "math"
  13. ;       subdirectory for the user-callable version of the 
  14. ;       incomplete beta function.
  15. ;-
  16.  
  17. function betacf, a, b, x
  18.   ;Continued fractions.
  19.   lc = a + b
  20.   ln = a - 1.0
  21.   lq = a + 1.0
  22.   max = 100
  23.   ja = 1.0 & jn = 1.0 
  24.   ka = 1.0 - lc * x / lq 
  25.   kn = 1.0
  26.   for i = 1, max do begin
  27.     ep  = i + 0.0
  28.     tm  = ep + ep
  29.     d   = ep * (b - ep) * x / ((ln + tm) * (a + tm))
  30.     jq  = ja + d*jn
  31.     kq  = ka + d * kn
  32.     d   = -(a + ep) * (lc + ep) * x / ((lq + tm) * (a + tm))
  33.     jq1 = jq + d * ja
  34.     kq1 = kq + d * ka
  35.     prev= ja
  36.     jn  = jq / kq1
  37.     kn  = kq / kq1
  38.     ja  = jq1 / kq1
  39.     ka  = 1.0
  40.     if(abs(ja - prev) lt 3.0e-7 * abs(ja)) then return, ja
  41.   endfor
  42. end
  43.  
  44. function ibeta_pdf, x, a, b
  45.  
  46.   on_error, 2  ;Return to caller if error occurs.
  47.  
  48.   if x lt 0. or x gt 1. then message, $
  49.     'x must be in the range: [0.0, 1.0]'
  50.  
  51.   ;gab = gamma(a) * gamma(b) 
  52.   ;gamma(a+b)/gab * exp( a*alog(x) + b*alog(1.0-x))
  53.  
  54.   if(x ne 0 and x ne 1 ) then temp = $
  55.     exp(lngamma(a+b)-lngamma(a)-lngamma(b)+a*alog(x)+b*alog(1.0-x)) $
  56.   else temp = 0.0
  57.  
  58.   if(x lt (a+1.0)/(a+b+2.0)) then return, temp * betacf(a, b, x)/a $
  59.     else return, (1.0 - temp * betacf(b, a, 1.0-x)/b)
  60. end
  61.  
  62.